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Abstract. Explicit expressions for arrival times of particles moving in a one- 
dimensional Zero-Range Process (ZRP) are computed. Particles are fed into the 
ZRP from an injection site and can also evaporate from anywhere in the interior 
of the ZRP. Two dynamics are considered; bulk dynamics, where particle hopping 
and decay is proportional to the numqber of particles at each site, and surface 
dynamics, where only the top particle at each site can hop or evaporate. We 
find exact solutions in the bulk dynamics case and for a single-site ZRP obeying 
surface dynamics. For a multisite ZRP obeying surface dynamics, we compare 
simulations with approximations obtained from the steady-state limit, where 
mean interarrival times for both models are equivalent. Our results highlight the 
competition between injection and evaporation on the arrival times of particles 
to an absorbing site. 



PACS numbers: 05.60.-k,87.16.Ac,05.10.Ln 



1. Introduction 

The Zero-Range Process (ZRP) is a stochastic model to describe the dynamics of far 
from equilibrium, interacting particles hopping between lattice sites [HO [3]. The ZRP 
has been used in many applications as a paradigm for transport processes, including 
traffic flows, shaken granular gases, network dynamics, phase separation, and particle 
condensation and clustering j4] . Mathematical interest also arises from the fact that a 
simple connection can be made between the ZRP and the totally asymmetric exclusion 
process (TASEP) pQ and that in certain cases - particularly for conserved systems - 
exact factorisable steady-state solutions can be derived [5]. 

In this paper, we compute the multiple passage times of particles obeying ZRP 
dynamics to reach a final absorbing site. We treat a nonconserved system where 
particles arc injected at the origin and evaporate as they drift right towards the end 
site of the lattice, as shown in Fig. [1] This type of dynamics may be applied to many 
specific micro-biological systems. For example, molecular motors may attach at one 
end of a microtubule, but desorb while traversing it. The distribution of arrival times 
of the motors will depend on their speed, and injection and desorption rates. Other 
examples include virus entry and transport to the nucleus, where the viral cargo is 
transported by molecular motors while being subject to degradation [6], and sperm 
entry into egg cells, where the first sperm to penetrate all layers of the cell triggers 
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Figure 1. Two realisations of a Zero-Range Process. In (a), any one of the red 
particles in the bulk can hop to the right or decay, while in (b), only the surface 
particle (in red) in each pile can hop or decay, the underlying particles being 
protected by the top particle. 



a block for subsequent ones [7J. In all applications there is a flux of particles, or 
"immigration," into the first site as well as particle annihilation at every location 
along the microtubule or layer within the egg cell. 

Unlike the first passage time (FPT) problem of a single, conserved particle 
undergoing a simple random walk, the first passage time of a nonconserved, 
multiparticle problem cannot be solved by analysing steady-states. In problems with 
injection and decay, the dynamics of particles reaching a specific "absorbing" site 
are complicated by subsequent particle arrivals, as well as arrival times conditioned 
on particles reaching the absorbing site. Since we will be concerned with a specific 
initial configuration, and wish to understand how the ZRP first reaches another 
configuration, we must find time-dependent solutions for the dynamics of the ZRP. 
Nonetheless, despite the nonconserved nature of the problem, steady-state solutions 
can still sometime provide useful approximations for FPTs of the ZRP in certain 
limits. 

We present exact solutions for arrival times in a finite-sized ZRP obeying two 
specific dynamical rules. In the first case, which we denote as "bulk dynamics" and 
which is illustrated in Fig. 1(a), all particles at a site are equally likely to hop to the 
next site. In the second "surface dynamics" case, depicted in Figure 1(b), only one 
particle can hop to its neighboring site. These two cases are limits of the ZRP and 
may serve as a model system for many physical systems. 

We begin our analysis in the following section with bulk dynamics. In this 
case, the particles are independent and we find exact analytic expressions for the 
distributions of the passage times for the k th particle to arrive at the absorbing site 
N + 1. For mathematical completeness, we present two derivations of the solution. 
The first involves explicit enumeration of the number of particles injected, evaporated, 
and having reached the absorbing site by time t. The second involves writing down 
a Master equation, which is solved using generating functions and the method of 
characteristics. All results for the bulk case are exact. In the third section, we define 
a ZRP with surface dynamics, where only one particle at each site, if it exists, is 
allowed to hop or desorb. In this case, we can only find exact solutions for a single 
site ZRP. For a multisite ZRP, we find particle arrival times in certain limits using a 
steady-state approximation, and compare them with results derived from Monte Carlo 
simulations. 
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2. Zero-range Model with Bulk Dynamics 



In our bulk-dynamic ZRP, beginning at time t = 0, particles are injected into the 
first site of an empty lattice whose positions are denoted by {1, 2, . . . , iV + 1}. The 
injection occurs as a Poisson process with rate a. Each particle can then hop one site 
to the right with rate p or evaporate with rate /i, independently of others. The forward 
hopping and evaporation processes continue until the particle reaches site N + 1 or is 
desorbed from the lattice. We wish to calculate the distribution of times for the k th 
particle to reach site N + 1. 



2.1. Explicit enumeration of particle fates 

Denoting by Tk the time at which the hnal absorbing site N + 1 is reached for the 
fc th time, we consider the accumulated number of hits H{t) by time t defined by 
P{H(t) =k) = P(T k <t< Tfe+i). The primary result of this section is that H(t) is 
Poisson distributed, with rate parameter 

n , n 



A(t) = aM-H*-^U 
w \n + p) V fi + pJ 



N-l . I , m 

n + p ^Kp + pj z^ o \ p + II ^ o n m \ y 

From this result we can find the survival probability Sk(t) that the final site has been 
hit by k — 1 or fewer particles, 

S k (t) = J2 P{H{t) = 3) = e~™ J2 • (2) 

3=0 j=0 ] ' 

Similarly, the distribution of 7* can be represented as the sum 

P{T { <*) = £ P(H(t) = 3) = e- A(t) £ ^f- (3) 

To derive the distribution of H(t), we begin by noting that we may break up the event 
{H(t) = k} according to how many particles n were injected before time t, each with 
common probability q(t) of reaching site N+ 1 by time t (see Fig. 2). The probability 
that exactly k of those n particles reach site N + 1 is Binomially distributed with 
parameter q(t), so that 



a ^ {f , +p)t ^f p Y^f n ^((/i + p)*)' ((p+pW 



n=k 



q(t) k (at) k „_ at g (a(l-g(t))tr- k 



-e 



k\ ^ (n-k)\ 

n=k 



= (gg^ e -a,(t) t- (4) 
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Figure 2. Space-time plot of ZRP bulk dynamics. In this realisation, 4 particles 
are injected by time t, and are the only ones capable of contributing to H(t). 
Trajectories that intersect the vertical line have arrived at site N + l within time 
t. Trajectories that intersect with the horizontal line at time t are those that 
failed to reach the absorbing site N + 1 by time t. 



This implies that H(t) is Poisson distributed with parameter 

A(f) = aq(t)t. (5) 

Deriving the probability q(t) is therefore our task. The particle injected at time r is 
characterised by its decay time £ r , which is exponentially distributed with mean 
and by its arrival time X T at the target site iV+1 excluding the possibility of decay. As 
this latter random variable is the sum of exponentials, it is Gamma(7V, fc)-distributed 
[9]. The probability q(t) of reaching site N + 1 is then given by the probability that 
the arrival time precedes both the chosen time limit t or the time of decay, averaged 
over the possible injection times r. Symbolically, 

1 r* 

«(*) = 7 / i P ( X - <t-T<(r)+ P(X T <(r<t- r)] dr. (6) 

* Jo 

Since X T and ( T are independent, the first probability P(X T < t — r < £ T ) = P(X T < 
t — t)P(( t > t — t). And since £ T is exponentially distributed with mean /i, and X T 
is Gamma-distributed, the first probability in the integrand of Eq.[6]is simply 



P(X T <t-T<Cr) = F r (t-r;N,p)[l- F Ct (t - r)] 

= F r (t-r;N,p)e-^- T \ (7) 

where Fr(s; N, /3) is the Gamma distribution function. The second probability, and 
many of the computations to follow, rely on the following equivalent representations 
of this function, 

gM. (8) 
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Here 7(2, w) = f Q t z ~ 1 e~ t dt is the lower incomplete Gamma function, and the 
right hand equality with the Erlang distribution holds because N is an integer [5]. 
Expression [5] leads to the following useful identities for F? : 



F r (u; N, /3)du 



1-e 



N-l 

7,T, 



-pu 



N-l 

E 

1=0 



N-l 



du = s 



E I /> 

^FrM + 1,/?) 



ufe'^du 



(9) 



and 



f e-^ F r (u; N, (3)du =-{l- e^ s ) - -J— ]T (-£-) V r (s; I + l,r? + /3). (10) 
Returning to the second probability in Eq.[6j Eq.[T0] yields 



P{X T <( T <t-r 



Jo 



^e^ s F r (s;iV,p)ds 



AT-l 



1 - e 



E(^T-) ^(t-r^ + l^ + p), (11) 



so that upon substituting Eqs. [71 and ITTI into Eq.[Sl we obtain 



1 r* 



Af-l 



/' 



- EH?") ^(t-r^+l.p + M) 



dr. 



Evaluating the integral termwise using Eqs.1^1 and 1101 we find 



1 N_1 



(xt + Fr{t;£+l,fJ,+p) 



F r (t;m + l,fi + p) 



m=0 



We can now expand the distribution functions in terms of a finite sum. Performing 
some algebraic simplifications yields the closed form representation 



N-l 



n + p 



H + p 



E 



(it 



((p + pW ((jM+ P )ty 



(12) 



m=0 r i=0 

Eq.[T2]may be further simplified by extracting the first two terms in the brackets from 
the sum and using the identity 



v^ fl , a(l-a n + 1 )-(n+l)(l-a)a n+1 
= 



1=1 



(i- a y 



(13) 
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Figure 3. Plots of the distribution P(H(t) = k) of the number of hits H(t) that 
have occured up to time t. Parameters used were p = 1, N = 10, fi = 0.2, and 
a = 3. The distribution is plotted for times t = 5, 10, 20, 30, 40, 50. Monte Carlo 
simulation (300000 runs) was used to verify the results for t = 30. 



obtained by differentiating the well-known identity ^?_ a £ = (1 — a n+1 )/ (1 — a). We 
finally end up with 



V/i + p/ V t(n + P )J t(p + P ) 



QW = 



x e 



h y i l+p) ^\p+»h il ml r 

which yields Eq.Q] In Fig. [3J and in all subsequent plots, we nondimensionalise all 
rates in terms of p and times in terms of p" 1 . The distribution P(H(t) = k) plotted in 
Fig. [3] shows that the probability of more arrivals increases with time after the start of 
injection increases. We compared and verified our results at time t — 30 with a Monte 
Carlo simulation of the bulk dynamics using the Bortz-Kalos-Lebowitz algorithm [12] . 
We performed 300000 runs, each with parameters t = 30, N = 10, p = 1, fi = 0.2 and 
a = 3. 

From the survival probability Sk(t) found from Eq. [21 all moments a of the fc th 
passage times can be computed, 

(2*) = -y f-^d*. (15) 
The mean (er = 1) first (fc = 1) passage time to the absorbing site can be found from 

p \ N ( N 



Si(t) = P(H(t) = 0) =exp 



^ + p/ V + P 



- 1 



M+P 



E(— )E — E 



(16) 



f=0 n m=0 r i=0 

and Eq.[l5] An explicit expression can be found for a single-site ZRP (N = 1): 
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(Ti) = — — — j(x,x), N = l, (17) 
fi + p 



where 

ap 



(18) 



Upon approximating the lower incomplete gamma function j(x, x) = T(x) — T(x,x) 
in the x — > limit, we find 



(Ti) = 



1 

H + p 



1 

x 



l + 0(x) 



N = 1. 



(19) 



In the x — > oo limit, we apply the method of steepest descents |16j to the integral 
definition of the T-function to find 



1 



H + p 



l + - + (D(x- 2 ) 



N = 1. 



(20) 



For N > 1, the integral in Eq.[l5]does not have a simple representation in non-integral 
form, nor can the mean k th passage times be calculated explicitly in the N = 1 case 
for k > 1. However we can find some asymptotic results for large k values in the "fast 
dynamics" limit. If we denote by r the characteristic time r = N/(fi+p) for a particle 
just injected to reach the final site, and consider times t such that t r Eq.Q]may 
be written as 

\(t) = a eS (t-T)+0(e- t ^), (21) 

where a e s = ap N (/j, + p)~ N is an effective injection rate from the perspective of the 
final site that takes into account decay. Because {> t holds for all but a negligible 
part of the range of the integral in EqfT5l we have 



(T k ) = — +r + 0(r 2 ). (22) 

Ctett 

The assumption on t translates onto a condition on k, so that Eq.[22] remains valid as 
long as (Tfc) » r, or k > a gr. 

In Fig. [H we plot the interarrival times (T&) — (Ife_i) as a function of k, fi, and 
ck. Fig. H{a) shows that larger desorption rates permit the system to reach steady- 
state faster, where interarrival times approach the limit l/a e ff, for smaller values of 
k. In Fig. B|b), we see that the mean interarrival times, including the mean first 
passage time, increase exponentially for large fi. This result is different from that of 
the problem where only a single particle is injected, in which the conditional mean first 
passage time of that single particle decreases has desorption fx increases. Because the 
single particle problem needs to be conditioned on arrival, only very fast trajectories 
can survive the desorption process, leading to mean first passage times that decrease 
rapidly with increasing /i. Finally, Fig. |4jc) plots the mean interarrival times as a 
function of injection rate a. 

In the next section we rederive the same results above by solving the corresponding 
Master equation using generating functions. Using this approach, we can not only 
recover the mean passage times to site N + 1, but the full particle occupation 
distribution function P(ni, . . . , njv+i, t). 
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Figure 4. Interarrival times (Tj.) — (Tfc_i) for the ZRP obeying bulk dynamics 
(with p = 1). (a) Interarrival times as a function of arrival k various values of 
the desorption rate fi for fixed a = 3. (b) Arrival times of the first, second, 
fourth, and sixth particles as a function of fj, for fixed a = 4. (c) First, second, 
fourth, and sixth arrival times as a function of the injection rate a for fi = 0.01. 
All plots were evaluated using N = 25. Results were verified using Monte-Carlo 
simulations (asterisks). 



2.2. Solving the bulk dynamics via generating functions 

In this section, for mathematical completeness, we rederive the survival probability 
using generating function methods applied to the Master equation for describing the 
probability P({n e }, t) of having {n e } particles on each of the 1 < I < N + 1 sites: 

P({n e },t)= -a[P({n^},t)-P(ni-l,...,niv + i,t)(l-5 ni ,o)] (23) 

N N 

-(»+p)J2 n J p (W}> i ) + ^E^' + ypfai ■ ■ ■ > n J + !» • • • > *) 

i=i i=l 

N 

+ P^2( n 3 + l ) p i n i,-- ■ ,nj + l,n j+ i - 1 ■ • ■ ,t)(l - S rij+u0 ). 
i=i 

The survival probability Si (t) is defined by the probability of having no particles in 
the absorbing site 



Si(t) = E P( n i>- ■ ■ > n N+i =0,i). 



(24) 



til ,"2---nw 



After setting tin+i = in Eg. 1231 we consider the dynamics of the constrained 
generating function defined as 



Go(zi, ■ ■ ■ , z N , t) 



roi,---n?v=0 



z^ N P{n u ---n N ,n N+1 =Q,t). 



(25) 
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Upon multiplying Eq.[23] by z™ 1 • • • z^ N and summing over all possible values of rij, 
1 < j < TV, we find a first order partial differential equation for Gq{z\, . . . , Zjf, t) which 
we can solve using the method of characteristics [T7] . We find that Go obeys 

= a( Zl - 1)G (26) 

at 

along the trajectories defined by 
dzi(t) 



dt 



= (fi + p)zi(t) - /J,- pz 2 (t), 



dzj(t) 

= {fJ- + P)ZjW- fJ--pZj+i{t), (27) 

dz N (t) 

— — = {fi+p)z N {t) - /i. 

The initial condition P(«i, • • • , n?v, Jijv+i = 0,t = 0) = S nit o ■ ■ ■ 5 nNt o gives 
Gq(zo, ■ ■ ■ ZN,t = 0) = f. Equations [27] can be written in the form 

^f=MZ(t)- M I, (28) 

where Z(t) = (zi(t), . . . , zjv(£)) t is the vector of trajectories, I is the N x N identity 
matrix, and M is a tridiagonal matrix with elements m,jj = \i + p, mjj+i = —p, and 
rriij = otherwise. Upon defining the Laplace transform Z(s) = Z(t)e~ st dt and 
the initial values Z(t = 0) = {z\(t = 0), ■ ■ • , z^it — 0)) T , Eqs.[25]can be written in the 
form 

sZ(s) = MZ(s)-^I + Z(* = 0), (29) 
s 

and solved explicitly by first inverting si — M and then calculating the inverse Laplace 
transform of Z(s). After performing the algebra, the solution to Eqs.^H] can be 
expressed as 

(Zj - Rif-j) = ]T (W* = °) ~ RN-j-k)^ L e^ t 1 (30) 

fc=0 



fe+i 

P 



where 

iZfc = l-[-*M . CM.) 

Upon using Zi(t) from Eq. [30] m Eq. [251 we find Go as a function of the Zj values 
implicitly expressed through the starting positions zj(t = 0) of the trajectorie: 



G (zi, ■■ ■ ,z N ,t) = exp 



J P \ N (z k +i(t = 0) - R N - k -i)p k 

- at {jT-p) - a g (p-Tpy^ki (32) 

x 7 [fe + l,-( M + ^] 
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We are thus left with explicitly determining zj (t = 0) as a function of the independent 
variables Zj. We do this by inverting Eq.1301 

N ~ 3 (r>t) k 
(*,(* - 0) - R,N-j) = £ (z 3+k Rn+jWJLe-^)*. (33) 



Using this result in Eq. [33] we find 



C?o(zi, • • • ,ZN,t) — exp 



-at 



P 



k=0 



N 



N-l N-k-1 



(»+ P )tY V pj+ " t3 (34) 



— ' (u + ») fc+1 A:!j! 
fc=0 j=o ^ ; J 



P + M 

x(z j+fc+1 - iijv-j-fc-x) 7[ fe + 1) -(M + P)* 



Finally, since the survival probability is obtained by imposing = 1 for all ^, we 
obtain 



S±(t) = exp 



-at 



P 



p + m 



2V 



ap e w^w 1 ^ ^ ir(ji+p) 
T Z^ Z^ 



(m+p) 



j=0 fcO 



7 [i + l,-(/x+p)t] 



(35) 



which is equivalent to Eq. 1161 We can now successively determine the dynamics of the 
probability distribution function conditioned on the absorbing site containing a finite 
number of particles n^r+i > 1. For rijv+i = lj the corresponding generating function 
Gi(z\, • • • , zn, t) can be obtained from the Master equation for P{n\, • • • , n^+i — 1, t): 



Gi(zi, • • • ,ZAr,t) = 2J ^■■■z^ N P(ni,---n N ,n N+ i = l,t). 



(36) 



ni,"-njv=0 



Upon summing Eq. 1231 over all values of rij we find that the dynamics of G\ is given 
by 



dGi dG 
a(j\(z\ — 1) + 



dt 



dzi 



(37) 



where Go is the generating function associated with the adsorbing site having no 
particles, n N+ i = 0, and where the evolution of the trajectories Zj(t) are unchanged 
from those described by Eqs.[57J The solution to Eq.[37]can be expressed in the form 



Gx(zi, ■ ■ ■ ,z N ,t) = X(t) G (zi, ■ ■ ■ ,z N ,t), 



(38) 



where X(t) obeys 



dA(t) _ p dG 
dt Go dzi 



(39) 
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The solution for X(t) turns out to be precisely that given in Eq.[T] Similarly, it can be 
found that the generating function with the constraint njy+i — j is given by 



' z N N P( n l, ■ ■ ■ n N+l = j,t) 

,z N ,t), (40) 

The survival probability Sk(t) is given by Sk(t) — Sj=o ^j0-> " ' j 1) an< ^ moments of 
the fc th arrival times, found previously, can also be obtained using Eq. 1151 In addition, 
an advantage of the generating function approach is that the particle occupations 
can also be determined. For example, the mean occupation at site i, conditioned on 
exactly j particles having entered site N + 1 is given by 



Gj( Zl> ---,z N ,t) = Yl z ^ 

ni,---njsf—0 

\(ty 



31 



-G (zi, 



(ne(t\n N+1 = j)) = ^ niP(m, ■ ■ ■ ,n N+ i = j,t) 

ni,---,njy— 

X(ty dG (l,---, z t , t) 



dzi 



'ap 



E 

k=0 



-k-l 



(fi + p) k + 1 M(£ - k - 1)! 



(41) 



Upon summing Eq. 1411 over all j, we find that the unconditioned mean occupation 
(rii(t)} is given by 



OO 

(n t (t))= Y n e P(n lr -- ,n N+1 ,t) (42) 
ni,---,njv+i=0 

iTo + k- 1)! ' 

Two limits are of interest: the mean occupation conditioned on no particles hitting 
site N + 1 , which is given by 



(n e (t\n N+1 = 0)) = e- A « (n e (t)), (43) 

and the average occupation of site N+ 1, regardless of the occupation state of all other 
sites, which is simply (n,N+i(t)} = X(t). Thus, in the long time limit, the occupation 
of the final site N + 1 will scale as 



(n w+1 (i))~crf ^- , (44) 

indicating that at site N+ 1 particles accumulate linearly at a rate that is proportional 
to the injection rate a attenuated by the evaporation probability for each of the N 
intervening sites. 
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Figure 5. Time dependence of the mean site occupancies. Both panels display 
exact values (solid lines) from Eq. 1431 and simulation (asterisks) for parameter 
values p = 1, N = 5, and fj, = 0.2. Each curve and approximating points 
correspond to mean occupations at different sites, with earlier site having higher 
occupations. 



In Figure [SJ we have plotted the mean occupations derived from Eq.03]for N = 5 
and fi = 0.2. All mean occupancies are seen to reach steady state by t ~ 10 and, for 
all times, mean occupancy is a monotonically decreasing function of site index due 
to decay. In Fig. 03a), a = 1 < p + /x, implying that particles are cleared out faster 
than they are injected, resulting in (rii) approaching a value less than one. In (b), 
the a — 1.5 > p + fi, and (m) (and (712) asymptotes to values greater than one. Our 
results are verified with Monte Carlo simulations. 

Finally, the full distribution for P(n\, • • • , n^+i = j, t) can be found by using the 
Cauchy integral [18] over Eq.[40l 



P(n u . . . , n N+1 = j, t) = — f ^+rd 2l • • • dzjy, (45) 

Z7rt Jc z 1 , • • • z N 

where the integral is closed along a path encircling the origin. Evaluating the residues, 
the above integral can be expressed as 



P( ni ,...,n N+l =j,t) = ^-l[( — ) Go( Zl ,---,z N ,t)\ ze=0 , (46) 

1=1 ^ Ze ' 

which can be calculated explicitly to yield 



P(m, . . . ,n N+ i = j,t) 



X(tyG (0,---,0,t) 



N 

n 



E 

k=0 



p l - L t 



lji-l-fc 



(p + At) fc+1 fc! 



T T f 7[fc + l,-0i + p)t] 



(47) 
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3. Zero-range Model with Surface Dynamics 

Surface dynamics differs from bulk dynamics in that only the top particle at a given 
site is able to hop or decay, while the ones below remain inert. The dependencies 
thus introduced between the particles renders the A:th hitting time, and even the 
survival probability Sfc(i) of the final site, too difficult to derive in the general N 
case. In particular, we cannot use the strategy employed for bulk dynamics because it 
relied essentially on particles injected before time t having independent probabilities 
of reaching site N + 1 by time t. In the case of surface dynamics, only the top particle 
in a pile attempts to move to a neighboring pile. The difficulty arises in keeping track 
of which sites are empty and which ones contain at least one particle. 

Beginning with the approach of the previous sub-section, we first consider the 
Master equation for the distribution of the site occupancies obeying surface dynamics 

P({n e }, t) = - a[P({n e ], t) - P(m - 1, ... , n N+1 ,t)(l - S nifi )} (48) 

N 

- n ]T[(1 - *„„o)P({n/}, t) - P(m, . . . , nj + 1, • • • , t)] 

N 

-pJ2(l-S nji0 )P({n e },t) 
i=i 

N 

~ P^-i 1 ~ s n j+1 , )P{ni, ■ ■ ■ ,Uj + l,n j+1 -l,...,t). 

3=1 

We now introduce the marginal probability 



Pi(rii,t)= P(n u ...,n N+ i,t), (49) 

where the sum is taken over all rij for all sites 1 < j < N + 1, except site j — i. Eq.l4l?l 
represents the probability that site i has ni particles regardless of the occupation of 
all other sites. Similarly, the joint probability for sites i — 1 and i is defined as 



Pi-i t i(ni-i,iii,t) = ^ P(ni,...,n N+ i,t). (50) 

» 3? !i-i,i=0 

Upon summing Eq.|48] over all values of rij^i, we find the time evolution for 
the marginal probability Pi(ni,t) as a function of the two-site probabilities 

Pi-i,i(ni-i,rii,t): 



Pi(rii,t) =p ^ (1 - 5„ i> o)P(n i _i,n i - l,t) -p ^ ^( n i-i> n i>*) + ( 51 ) 

n-i — i — 1 rii-i— 1 

+ (/j, + p)P l {n l + 1, t) - ( M + p)P(m, t)(l - ^,0), K i < N + 1 

(52) 

Continuing in this way, the equations for the marginal occupation probabilities form 
a hierarchy which is completed by the equation for the injection site i = 1: 
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A(m,*)= -a[P 1 (ni,t)-P 1 (ni-l,t)(l-5„ 1)0 )] (53) 
- ( A t + p)[Pi(m,i)(l - 5 ni , ) - Px{nx + l,t)]. 

Note that the dynamics for site i = 1 is completely decoupled from that of the other 
sites so that the marginal occupation distribution of the first site can be solved directly. 
We now consider two cases where analytic results can be found. 

3.1. Single- site ZRP densities and mean first passage times 

Since Eq.[53] is decoupled from the hierarchy, we can solve it by taking its Laplace 
transform and using the initial condition P(m,0) = <5 ni .o to find 



^ , s + a ~ , „ 1 

Pi m = 1,8) = — — Pi Til - o,s) — 

l i +P M + V 

t S -\- CX\ ~ o 

Pi(rn + l,s )= 1 + — Pi(ni,s) — Pi ni - l,s). 

V M + P/ M + P 



(54) 



The solution can be expressed in the form 

1 - zx(s) 



zi{s) n \ (55) 



Px(nx,s) 
where 

ZiO) = yj^pj {{s + a + fjL+p)- v / (s + a + (i+p) 2 -4a(/i+p)) ■ (56) 
Upon inverting, we find 

zx(t) = - __ V h(2^a^ + p)t), (57) 

where -Zi(£) is the first-order modified Bessel Function of the first kind. From Eq.l55l 
and using the fact that the inverse Laplace transform of a product is a convolution in 
time, we can iteratively construct Pi(n\,t) starting from Pi(0,i) 

Pi(ni = 0,t) = 1 - [ zx(t')dt', (58) 



where z\ (t) is given by Eq.[57] and 



Px(nx,t)= f Px(nx - l,t')zx(t')dt' 
Jo 



(59) 



The integrals in Eq.[SH]and[5S]do not have simple closed forms. However, the functions 
Pi(ni,£) can also be obtained from differentiation using the relation 

Px(nx + l,t) =( i-5„ 1<0 + _2_) Pi ( ni ,t) (60) 

| P 1 {n 1 ,t)-aP 1 {n 1 -l,t)(l-6 ni , ) 
M + P 
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For instance, we may recursively write 

= l,t) = (l - [ t z l (f)dA - (61) 

In the case of N ~ 1 we can also solve for the first (k = 1) passage times by observing 
that the equation for the two-site distribution function, conditioned on 712 = 0, is also 
decoupled from the hierarchy: 

P(m, 0, t) = a[P(m - 1, 0, t)(l - <5„ 1)0 ) - P(ni, 0, *)] + M^K + 1, 0, t) 
-( f i+p)P(n 1 ,0,t))(l-S nu0 ), 

where we have dropped the subscripts on the two-site distribution function so that 
P(n\,ri2,t) = Pi,2(ni,ri2,t). Using Laplace transforms, we find 

P(m,o,s)= yi(s)ni (62) 

s + a - nyi(s) 

where 

a + p + p + s - y/(a + p + p + s) 2 - Aap 

yi(s) = ^ . (63) 

From this solution of P(ni,ri2 — 0,s), we can obtain the Laplace transformed 
probability that site i = 2 has not been hit by a particle £>i(s) = X^=o P( n ^i n 2 = 
0, s). Thus, in the N = 1 case, the mean first passage time is 



irp \ £, n nA a + P + P+ yj{a + p + p) 2 - lap 

(Ti) = 2^ p ("i'°' s = °) = 2 ■ ^ ' 

Note that in the case of p — 0, this result simplifies to (Ti) = a -1 In surface 

dynamics without desorption, the first passage time is determined by the dynamics 
of the lead particle. Therefore, the mean first arrival time is simply the total time it 
takes for the leading particle to reach site N + 1 and is given by (Ti) = or 1 + Np~ x . 

3.2. Steady- state limit 

We have not been able to find closed-form solutions of the surface dynamics ZRP for 
general N and nonzero desorption rate p > 0. However, Eq. [53] can be solved in the 
steady state limit by using the ansatz P\{n\) = (1 — z\) z" 1 . The equation supports a 
solution when z\ = a(p + p) _1 , implying 



P l (m)= 1--=- — • (65) 

The above expression is correct only for a < p+p, so that < Pi(ni) < 1. Physically 
this condition is simply a statement of that if injection is too fast, the occupations 
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continues to build without bound. Steady-state occupations arise only if the injection 
rate a is small enough such that hopping p and evaporation /i can keep up. 

In order to solve Eq.(2U we need to a closure relation for the two-site probability 
distribution Pi— ii{TH—i, rij). As shown in [5], the two-site probability distribution 
can be factorised in the steady-state limit and expressed as Pi_i^(rij_i, rij) = 
Pi-i(ni-i) Pi(ni). If we impose that each Pj(rij) has a power law dependence in 
nj, similar to what done for m, it is easy to verify that the steady-state marginal 
probabilities are solved by 

The resulting steady-state mean occupation at each site are 



= 7 — —r: —, a<^+p. (67) 

[ji + py - ap 1 

From our steady-state results for N > 1 , we can find an approximation to the passage 
times by a mean-field argument in which the probability of site N + 1 surviving 
up to k — 1 particles hitting it obeys S'fc(i) = —J(t)Sk(t). The particle current 
J(t) = p^] j P(nN\Tk > t) is conditioned on fewer than k particles having arrived 
at site 7V+1 by time t. Since neither P(njv|Tfc > t), nor the unconditional distribution 
P(riN,t) are available, we must approximate J(t) with its steady-state, "mean-field" 
(single site marginal distribution) value through Eq.l6"6l 

oo 

J Rip P(riN, t — > oo) 

"jv = 1 

N — 1 \ / N-l 

ap \ I ctp 



JV 

ap 

In this approximation, J is independent of k and all interarrival times are 
approximately 

(T k ) (n-t) = / S k (t)dt « (69) 
Jo ap N 

As expected, this estimate is precisely that given by a e ff in Eq.[2T] for the bulk 
dynamics case and is accurate in the limit of a <C (/i + p) where the entry flux is 
slow compared to the internal dynamics and the first passage time is dominated by 
the contribution given by a" 1 . Fast internal dynamics allows the system to quickly 
reach steady-state, rendering the interarrival times equivalent for bulk and surface 
dynamics. 

Upon taking the limit of slow injection rate a — > in Eq.|B3]we find 



limm) = 0i±5>, (70) 

a-s-0 ap 

which is identical to the result in Eq.[5H] for N = k = 1. Figure |5] plots simulated 
interarrival times and compares them with those from bulk dynamics. For large k, 
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Figure 6. Interarrival times (Tj.) — (Tfc_i) for the surface dynamics ZRP as a 
function of k. Results from Monte Carlo simulations (open symbols) for both 
N = 1 and N = 10 are presented. Since fewer particle are mobile in surface 
dynamics, the arrival times are longer. The simulations match the analytic results 
found for the N = k = 1 (Eg , 1641 indicated by the asterisk) and large k (Eqs. l69l 
and !22ll limits. For comparison numerical results for the bulk dynamics case (filled 
symbols) are also plotted. All plots were derived using fi = 0.1 and a = 1. 



interarrival times for both bulk and surface dynamics approach the same value a~ s 
for each N. The only exact result for surface dynamics is that given by Eq. 1641 for 
N = k = 1, and is indicated by the asterisk at (Ti) w 2.05125. 

4. Summary and Conclusions 

In this paper, we have provided detailed and explicit calculations of first passage 
times of an N— site, one-dimensional Zero- Range Process. Both a Poissonian injection 
process at an injection site, and spontaneous desorption of all sites were included. We 
considered both bulk dynamic and surface dynamic rules as illustrated in Fig. [TJ 

For the ZRP obeying bulk dynamics, we computed the particle passage times 
using two methods. In the first method, we explicitly enumerated the random walks 
of each injected particle and evaluated their probability of reaching the final absorbing 
site within a time window. The probability that the absorbing site has not absorbed 
more than k particles by a certain time was constructed. The main results for the 
survival probabilities are given by Eqs.[2]and[I] with explicit expressions for the mean 
first passage time given by Eq.[17]and its subsequent asymptotic limits. 

We also derived the complete Master equation for the probability distribution for 
a ZRP obeying bulk dynamics, and solved its corresponding generating function using 
the method of characteristics. In addition to the k th passage time distribution, this 
yielded the mean conditional occupancies of each site given by Eq.UI] and the full 
probability distribution given by Eq.l47l 

Finally, for a single site (N = 1) ZRP obeying surface dynamics, we found exact 
results for the site density distribution (Eqs. [5S] and I59|) and the mean first passage 
times fEq. lM|) . Note that higher moments of the first passage time are readily obtained 
by evaluating higher derivatives of Eq.[52]at s = 0. For general N, only the steady- 
state particle currents and interarrival times could be found in closed form (Eq. l69[) . 
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